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ABSTRACT 

In  this  study  a  numerical  model  based  on  the  hydrostatic  Boussinesq  equations  is 
used  to  simulate  atmospheric  frontogenesis  driven  by  an  irrotational  non-divergent 
deformation  wind  field  The  equations  are  numerically  mtegrated  by  using  the  serru- 
Lagrangian  technique  associated  with  two  different  time  schemes:  explicit  and  senii- 
implicit  Both  schemes  produce  realistic  fronts  after  approximately  40  hours  of  model 
integration.  Thesemi-Lagrangian  semi-implicit  scheme  is  more  successful  in  handling 
the  sharp  gradients  associated  with  the  front.  Also,  the  semi-Lagrangian  semi-implicit 
equations  are  mtegrated  with  time  steps  as  long  as  3600  sec.  producing  solutions  with 
relatively  small  errors.  This  indicates  that  this  numerical  scheme  is  appropriate  for  use 
in  mesoscale  regional  models. 
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I   INTRODUCTION 

The  development  of  fast  computers  has  been  of  great  importance  for  the 
atmospheric  sciences.  With  them,  scientists  have  been  able  to  obtain  a  better 
understanding  of  the  physical  processes  in  the  atmosphere  and  operational  weather 
services  have  produced  more  accurate  forecasts.  But  as  pointed  out  by  Robert  ( 198 1 ), 
scientists  can  also  give  their  contribution  by  developing  numerical  algorithms  that  are 
computationally  less  expensive  and  at  the  same  time  accurate  This  interaction  between 
science  and  technology  will  allow  the  computational  effort  to  be  used  in  numerical 
models  that  will  have  higher  resolution  and  that  will  be  more  accurate,  as  well,  since 
more  realistic  and  complex  physical  processes  would  be  represented  in  these  models. 

One  of  the  atmospheric  processes  that  has  been  under  active  study  is  the  formation 
of  fronts.  In  a  study  of  non-geostrophic  frontogenesis  Will iams(  1972)  reproduced 
numerically  a  realistic  front,  demonstrating  the  importance  of  nonlinear  effects  in  this 
dynamical  process.  Since  frontogenesis  is  a  phenomenon  that  is  characterized  by  the 
development  of  sharp  gradients  of  the  physical  parameters  it  would  be  interesting  to 
develop  a  numerical  model  based  on  techniques  that  can  handle  adequately  the  scale 
collapse  problem.  Kuo  and  Williams  (1990)  studied  the  use  of  different  numerical 
techniques  in  the  solution  of  a  simple  scale  collapse  problem.  They  concluded  that  the 
semi-Lagrangian  method,  developed  by  Sawyer  ( 1963),  would  be  adequate  for  use  in 
numerical  models  where  strong  gradients  were  present. 

The  semi-Lagrangian  technique  has  been  used  mostly  in  synoptic  scale  numerical 
models.  Robert  ( 198 1)  used  this  technique  with  the  semi-implicit  scheme  ( Robert  et  al.. 


1972)  in  a  primitive  equation  model  and  obtained  perfectly  stable  solutions  using  time 
steps  25  times  larger  than  those  that  would  be  alJowed  by  the  Courant-Fnedrichs-Lewy 
(CFL)  stability  criterion.  Pudykiewicz  and  Staniforth  ( 19X4)  examined  the  use  of  the 
scheme  in  several  different  conditions  and  concluded  that  besides  its  good  stability 
properties,  it  is  also  accurate  and  flexible. 

The  objective  of  the  present  study  is  to  perform  experiments  using  the  semi- 
Lagrangian  technique  in  a  frontogenesis  model  and  verify  its  suitability  for  use  in  the 
representation  of  atmospheric  processes  where  sharp  gradients  are  present 

The  semi-Lagrangian  technique  is  introduced  in  the  next  chapter.  The  basic  model 
equations  are  developed  in  chapter  III  The  numerical  procedure  employed  for  solving 
the  equations  is  presented  in  chapter  IV.  In  chapter  V  the  experiments  and  results  are 
discussed,  followed  by  the  summary  and  conclusions  in  chapter  VI. 


II.  THE  SEMI-LAGRANGIAN  SCHEME 

This  model  uses  a  semi-l  agrangian  algorithm  similar  to  the  one  described  by 
Robert  (1981). 

The  three  time  level  serru-Lagrangian  technique  consists  of  evaluating  the  total 
derivative  of  a  general  dependent  variable  Q(F(f),t)  following  the  trajectory  of  the 
fluid  particle,  using  the  approximation 

dQ  =  <2(f(r+AQ,r+Ar)-Q(r(t-Af)v-At)  (2.1) 

dt  2Af 


where  Fit)  represents  the  position  vector  of  the  fluid  particle  at  time  t  For  sake  of 
simplicity,  constrain  the  problem  to  two  dimensions,  y  and  z.  In  this  case  the  position 
vector  will  be 

m  -  y(oAz(o*.  I22) 

The  location  of  the  particle  at  the  forecast  time  (t  +  A  t)is  chosen  to  be  coincident 
with  a  grid  point.  Let  the  displacements  of  the  particle  along  the  y  and  z  directions 


during  one  time  step  be  represented  by  a  and  b .  respectively    If  the  grid  point  position 
is  represented  by 

r(r  +  Af)  -  (y.,zk)  (W) 

the  approximation  (2. 1)  can  be  written  as 

dQ  m   Q(yj,zk,t  +  At)-Q(yj-2a9zk-2b,t- &t)  (2.4) 

~dt  2Af  ' 

The  displacements  a  and  b  can  be  calculated  iteratively  using 

an+1  -  At.v(yra\zk-b\t),  <2-5a) 


bn+l  =  Af.w(v.-an,zt-ftn,r),  <25b> 

where  the  superscript  n  represents  the  n  iteration  and  v  and  w  represent  the 
horizontal  and  vertical  velocities,  respectively.  The  initial  estimates  of  a  and  />in  the 
iterative  process  are  defined  by 


a0  =  At.v(yjyzk,t), 


(2.5c) 


b°  -  ±t.w(yj,zk,t) 


(2.5d) 


In  Pudykiewicz  and  Staniforth  ( 1984)  it  can  be  found  that  the  necessary  condition  for 
convergence  of  the  iterative  procedure  described  in  (2.5)  is 


Af .  max 


V 

dv 

ay 

» 

dv 

dz 

> 

dw 
dy 

» 

dw 

dz 

<\ 


(2.6) 


Furthermore.  Kuo  and  Williams  ( 1 990)  point  out  that  no  more  than  three  iterations  are 
necessary  for  convergence  to  be  attained.  Robert  ( 198 1)  compared  the  use  of  two  and 
four  iterations  in  his  primitive  equation  model  and  the  results  obtained  were  quite 
similar. 

In  order  to  determine  the  values  of  the  variable  Q  at  time  t-At  it  is  necessary  to 
use  an  interpolation  scheme.  In  this  model,  bicubic  spline  interpolations  (de  Boor.  1 962 ) 
were  used  with  the  explicit  formulations  of  the  semi-Lagrangian  scheme,  using  the 
algorithm  described  by  Marchuck  (1982)  In  the  semi-implicit  formulation  only  one- 
dimensional  cubic  splines  in  the  y  direction  are  necessary  to  solve  the  problem;  as  will 
be  shown  in  [V.C.2. 


III.   BASIC  EQUATIONS 

This  model  uses  the  hydrostatic  Boussinesq  equations,  which  assume 
incompressibility  of  the  atmosphere.  Friction  and  heating  are  neglected  It  is  also 
assumed  that  the  atmosphere  is  bounded  above  by  a  rigid  lid.  Periodic  boundary 
conditions  are  used  in  the  horizontal. 

The  Boussinesq  equations  can  be  written  as 

*L   .    -**♦*,  (3.1) 

dt  dx 


dt  dy 


<*  +  <*+<>»   .Q<  (33) 


dx     dy     dz 


**    -  0  ,  (3-4) 

dt 


(3.5) 


where 


K  -  —  ,  (3.6) 


e  -  n£-)  -e0,  (3.7) 


<P   '  CB0(-£-)\*z.  (3.8) 


This  model  assumes  that  the  Coriolis  parameter  /is  constant,  and  that  the  bottom 
topography  is  flat.  The  variations  in  the  x  direction,  along  the  front,  are  neglected, 
except  for  the  «p  field  and  the  basic  deformation  field. 

With  the  assumptions  above,  equation  (3.3)  reduces  to 

f*+-^=0  (3-9) 

dy     dz 


for  departures  from  the  deformation  wind  field. 


With  the  use  of  a  rigid  lid  assumption  the  pressure  is  not  known  at  any  level,  and 
this  will  require  some  manipulation  of  the  equations. 

Integrate  the  hydrostatic  equation  (3.5)  from  z  -  #to  an  arbitrary  level  z  and 
obtain 


<p  =  j&  dz  +  C.  0-W) 


0% 


The  constant  of  integration,  C,  can  be  determined  by  taking  the  vertical  average  of 
equation  (3. 10)  and  subtracting  from  equation  (3.10)  which  gives 


<p  =  -£-[Qdz-<-^-[Qdz>  +  «i», 


(3.11) 


where  <  >  denotes  the  vertical  average  operator 


<F>  =  IfFdz.  <312) 


The  constant  of  integration  has  been  eliminated,  but  now  an  expression  for  <f»>must 
be  found. 


Apply  the  vertical  average  operator  to  the  continuity  equation  (3.9)  and  use  the 
vertical  boundary  conditions 

w  =  0       at       z  =  0  (3  13a) 


and 


w  -  0       at       z  -  H  tf-13b* 


which  vields 


^   =  0.  (3.14) 

3y 


Next,  expand  the  total  derivative  in  equation  (3.2),  multiply  (3.9)  by  v  and  add 
to  (3.2)  to  obtain 

dv     dw     dwv  dip     j.  ,,  ,-v 

+ + JL-fu.  (3.15) 

dt      dy        dz  dy 

Take  the  vertical  average  of  (3.15),    using  the  vertical  boundary  conditions  (3.13), 
which  gives 


3<v>  |  a<w>  _  _d<s>_f<u>  (316) 

dt  dy  dy 


Differentiating  (3.16)  with  respect  to  y  and  using  equation  (3. 14)  yields 

&«$>         a2<w>    rd<U> 


f^-=-  .  (3.17) 


dy2  dy>  fy 


With  appropriate  boundary  conditions,  a  solution  for  <^>>can  be  found 

In  order  to  solve  equation  (3.1),  an  expression  for  — —    must  be  obtained. 

dx 

Differentiate  (3.11)  with  respect  to  x  and  use  the  assumption  that  6 is  not  a  function  of 
x  which  gives 

dy   m   6<cp>  (318) 

dx  dx 


An  expression  for  — —  must  be  obtained.  Take  the  domain  average  of  (3.2) 
cbc 


{<^>}  -  {<-^P>}-{</I/>},  (3.19) 

dt  dy 


where  {( )}  denotes  the  horizontal  average  operator 


-j!Fdy 


(F>  -  ±f  Fdv.  (32°) 
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Expand  the  first  term  of  (3.19)  and  use  the  vertical  boundary  conditions  (3. 13)  which 
yields 

{<dv>}  m    di<v>\  m    d<v>  (U1) 

dt  dt  dt 

Using  periodic  boundary  conditions  in  the  y  direction,  it  can  be  shown  that  the 
second  term  in  (3.19)  is  zero.  With  the  assumption  of  constant  f,  and  using  (3  21), 
equation  (3 .19)  becomes 

<3<v> 


dt 


-   -f[<u>).  (3.22) 


Likewise,  take  the  domain  average  of  (3. 1 )  and  use  the  vertical  and  horizontal  boundary 
conditions  which  yields 

di<U>)            a(<<p>)      _    ^  (3  23) 
= X — +/<V>  . 

dt  dx 


Equations  (3.22)  and  (3.23)  show  that  if  {<u>}  is  set  initially  equal  to  zero.  <v> 

will  remain  constant  in  time  and  so  will  3- —    Furthermore,  if  <v>  is  also  set  to 

zero  initially, —  will  be  identically  zero  for  all  time. 

dx 


11 


Ln  this  study  it  will  be  further  assumed  that  — — —  is  constant  in  the  cross  frontal 

dx 

direction  (y)y  which  allows  to  write 


di<(p>}         3<cp> 
dx  dx 


With  the  assumptions  above  the  following  relation  is  obtained 


(324) 


d(p        d<cp>    h  ~ 


(3.25) 


for  initial  conditions 

{<u>'t  =  <v>  =  0.  <3-26) 


The  dependent  variables  will  be  expressed  as  a  combination  of  a  basic  and  a 
disturbance  part  as  follows: 

u(x,y,z,t)  -  £/,(x)y)  +  «/0>,z,r),  <3  27> 


v(x,y,z,0  -  K(,(x,y)  +  v'(y,z,l),  <•««) 


w(y,z,t)  -w'(y,z,t),  °29) 
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6(y,z,r)  -  #(y,z9t), 


n  3m 


<p(*,y,z,0  -  *(jc,y)  +  9/(y>z.O. 


(3.31) 


where  the  quantities  with  primes  represent  the  disturbance  part  of  the  variables 

The  basic  wind  field  that  drives  the  frontogenesis  is  the  nondivergent  and 
irrotational  deformation  field  given  by 

Ud(*>y)  '  -  —  sinh(^x)sin(^y),  0-32) 


VAx>y)  -  -—  cosh(^Jc)cos(n;y).  <-V33) 


This  wind  deformation  field  was  originally  used  by  Stone  (1966)  in  his  study  of  quasi- 

geostrophic  frontogenesis. 

It  has  been  assumed  that  the  perturbations  are  independent  of  x     U  the 

expressions  for  u  and  v  (3.27  and  3.28  respectively)  are  substituted  into  the  basic 

equations,  —  and  —  will  be  functions  of  x.  In  this  study,  where  the  main  interest  is 
dt  dt 

in  examining  cross  front  variations,  the  equations  will  be  evaluated  at  x  -  0.  This  will 
make  the  mathematical  formulations  compatible  with  the  previous  assumptions  and  it 
is  expected  that  the  results  will  be  physically  consistent. 
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Therefore,  the  basic  wind  deformation  field  will  become 

UA0,y)  -  0  ,  (3.34) 


Vd(.0,y)  -  -*™cos\iy.  0.35) 


The  basic  geopotentia!  field  #  is  chosen  such  that  it  will  be  in  geostrophic  balance 
with  the  wind  deformation  field  as  follows: 

fUd-   ~,  (3.36) 

dy 


~fVd-  "^-  (337) 

ax 


Substitution  of  (3.36)  and  (3.37)  into  (3.27)  and  (3.28)  respectively,  and  use  of  the 
expressions  obtained  for  the  dependent  variables  in  equations  (3.1)-(3.4)  yields  the 
following  set: 


--/V,  (338) 

dt 
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dv  dip     f      ,      f/  A  dVd 


fU-(v+Vd)-^t  (3-39) 


dt  dy  a     dy 


^+^    -0,  (3.40) 

dy     dz 


^-0.  (3.41) 


where  the  prunes  were  dropped  from  the  disturbance  variables. 

da 
Finally,  an  easier  way  to  obtain  —  will  be  described.  The  dependent  variable  q> 

dy 

can  be  expressed  as 


<P  -  <Po+<P,> 


(3.42) 


where 


•.-/£*• 


»eo 
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Differentiate  (3.42)  with  respect  to  y,  to  obtain 


1*    .   ^0^8_f3&dz  (3.44) 

J      rhi 


dy      fy    Vn  ^y 


The  value  of  ^is  not  known,  by  virtue  of  the  rigid  lid  assumption.  In  practice  v 
will  be  predicted  using  equation  (3.44)  with  — -  assumed  initially  equal  to  zero.  After 

ay 

that,  the  v  field  will  be  adjusted  in  order  to  satisfy  <  v  >  -  0.  This  procedure  will  be 
equivalent  to  solving  equations  (3.11)  and  (3.17)  for  tp,  differentiating  the  values 
obtained  with  respect  to  y,  and  substituting  into  (3.39).  Equations  (3  38)-(3  41)  and 
(3  44)  give  a  complete  set  that  can  be  integrated  in  time  in  order  to  predict  the  values  of 
the  dependent  variables. 
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IV.    NUMERICAL  SOLUTION 


A       INITIAL  CONDITIONS 

The  initial  conditions  used  in  this  model  are  similar  to  the  ones  used  by  Williams 
et  al.(  1991)  in  their  study  of  effects  of  topography  on  fronts. 
The  initial  temperature  field  is  given  by 


dz  2 


(4.1) 


where 


c  -  A  cosCji^cos2 


V 

\i.> 


71 

~2 


zzzt 


(4.2) 


e  -  0, 


z>z. 


(4.3) 


ae 

In  the  expression  above  — -  is  a  constant,  A  is  the  amplitude  of  the  temperature 
disturbance  and  z,  is  the  height  of  the  upper  limit  where  the  temperature  disturbance  is 
present.  This  initial  temperature  field  will  confine  the  frontogenesis  to  the  lower  layers 
of  the  atmosphere,  which  will  make  the  results  more  realistic. 
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The  initial  u  component  field  is  chosen  to  be  in  thermal  wind  balance  with  the 
temperature  field  and  is  given  by 


(4.4) 


where 


A.  -  (z-z,)-—  sin 


7t 


Z±Zt 


(4.5) 


and  the  condition 


has  been  used. 


k  -  0, 


z>zt 


(4.6) 


u(y,zt,0)  -  0 


(4.7) 
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The  initial  ^and  w  fields  are  obtained  from  the  solution  of  the  quasi-geostrophic 
circulation  equation  (AVilliams,  1972) 


/2e0    dz    dy2  +  dz2       /2e0    ay    fy 


where 


(48) 


V  -    -  **  (4.9) 

dz 


W  -   ^4-  (4.10) 

ay 


The  vertical  boundary  conditions  for  the  solution  of  (4.8)  are  set  to 

♦  (z  -  0)  -  \|r(z  -  #)  -  0  (4-n> 


It  can  be  shown,  from  the  domain  averaging  of  equation  (4.2)  that  {<  w  > }  will 
be  zero  initially.  Also,  the  condition  (4.11)  will  allow  (3.25)  and  (3.26)  to  be  satisfied  in 
order  for  the  numerical  solution  to  be  consistent  with  the  assumptions  of  the  model 
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B       THE  SPATIAL  GRID 

The  numerical  calculations  are  performed  on  a  staggered  grid,  both  in  the 
horizontal  and  in  the  vertical  directions.  In  the  horizontal  the  variables  are  staggered 
following  the  scheme  B  described  in  Arakawa  and  Lamb  (1977).  This  scheme,  as  shown 
in  Haltmer  and  Williams  (1980),  allows  a  better  representation  of  geostrophic 
adjustment  processes.  The  variables  are  also  staggered  in  the  vertical.  Pielke  ( 1984) 
points  out  that  the  staggering  of  the  u  and  w  components  in  the  vertical  gives  better 
solutions  for  the  vertical  velocity  from  the  continuity  equation. 

The  horizontal  grid  uses  I  grid  points  separated  by  a  constant  grid  space  Ay.  The 
vertical  grid  has  K  grid  points  and  also  uses  a  constant  grid  spacing,  referred  to  as  Az. 
Figure  1  shows  the  arrangement  of  the  variables  on  the  grid. 

C.      FORMULATION  OF  THE  MODELS 

1.      Semi-Lagrangian  Explicit  Scheme 

In  the  semi-I  agrangian  explicit  formulation  of  the  present  model,  the  total 
derivatives  will  be  approximated  using  the  technique  described  in  chapter  II.  The 
remaining  terms  of  the  equations  are  evaluated  at  the  positions  of  the  fluid  particles  at 
time  level  /  Introducing  those  approximations  into  the  prognostic  equations  (3.38), 
(3.39)  and  (3.41)  the  following  set  is  obtained,  after  dropping  subscripts  /andy : 

u(y,z,t+At)-u(y-2«yz-2b,t-*t)  _/v(  _M>  (4,2) 

2A; 
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w  •  w  •   w     z(k  +  l) 


e,<p  •  •  e,q>  •  •  e,<p    z(k+*; 

u,v,Vd  u,v,Vd 


w  •   w  •   w  z(k) 


e,(p  •  •  e,9  •  •  e,<p    z(k-i) 

u,v,Vd  u,v,Vd 


w  •   w  •   w        z(k-l) 


y(j-D  y(j-i)  y(j)  y(j+i>  y(j+D 


Figure  1:  Staggered  spatial  grid 
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v(y,z,f+Ar)-v(y-2a,z-26,r-  AQ 
2Af 


-^-(y-a,z-b,t)-f.u(y-a,z-b,t) 

dy 


[v(y-a,z-b,t)  +  Vd(y-a)]  -^(y-a),  (4.13) 


d(yyz,t+Af)-Q(y-2ayz-2b,t-At)  _  Q  (4  14) 

2Af 


The  values  of  wand  <p  are  obtained  from  the  diagnostic  equations  (3.40)  and  (3.45), 
respectively. 

2.      Semi-Lagrangian  Semi-Implicit  Scheme 

The  semi-Lagrangian  semi-implicit  scheme  was  introduced  in  this  model 
following  the  development  described  by  Robert  et  al.(1985).  The  scheme  consists  of 
separating  the  temperature  into  a  basic  state,  dependent  only  on  the  vertical  coordinate 
z,  and  a  disturbance  part.  Next,  an  implicit  treatment  is  given  to  those  terms  related  to 
gravity  waves,  namely:  horizontal  pressure  gradient,  divergence  and  vertical  motion. 

To  obtain  the  system  of  equations  in  the  semi-implicit  formulation  take  the 
advective  terms  out  of  the  total  derivatives  in  equations  (3.38),  (3.39)  and  (3.44); 
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d*U             ZU  t  ,4K) 

-  +  W  -  /V,  14.15) 


dt  dz 


dt         dz  dy  dy 


dt  dz 


where  the  terms  with  subscript  H  represent  horizontal  components  of  the  total 
derivative. 

Replace    the   temperature  and  geopotential    fields   by   basic  state  and 
disturbance  parts  as  follows: 


e(y>z,t)  -  e*(Z)+e/(y(z,r))  <418> 


<p(y,z,«)  -  v'(z)  +  <fXy,z,t),  ,419> 


where  the  superscripts  (*)  represent  the  basic  state  parts  ol  the  variables. 

Now  introduce  (4. 18)  and  (4. 19)  into  equations  (3.5),  (3.40)  and  (4. 1 5)-(4. 1 7), 
and  take  the  time  average  of  those  terms  related  to  gravity  waves,  which  yields 
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dHu  du       A 

dt  dz 


dz         80 


where  (     )  denotes  the  time  average  operator. 


(4.20) 


dt         dz  dy  d    dy 


dy      dz 


d$  dQ'     -dd*        n  (4  21) 

dt  dz  dz 


^£1    -    &L  (4.24a) 

dz  6° 


(4.24b) 
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Define  the  explicit  terms  at  time  t  as 


r     =    vv—  -/V,  (4  25) 

1        & 


r,  =  w^  .  <4-27> 

3        az 


Usmg  (4.25)-(4.27)  in  (4.20M4.23),  the  following  set  of  equations  is  obtamed: 


V 

<* 


+ 


r,  -  0, 


(4.28) 


^  +  i*+,    -  0,  (4.29) 
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^  +  -^    =  0  ,  (4.30) 

By      dz 


dt  dz      3 


Next,  define  the  auxiliary  variables 

F+  -  F(y,z,*+Af),  (432a> 


F°  -F(y-a,z,f),  (432b) 


F-  =  F(y-2a,z,r-Af), 


where 


(4.32c) 


an+1  -  Ar.(v°(y-flrt,r)  +  F/(y-aB,r)).  (4.33) 

Note  that  the  use  of  the  time  average  of  the  vertical  velocity  w  eliminates  the 
need  of  considering  vertical  displacements  in  the  calculation  of  the  trajectories  of  the 
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fluid  particles.   Therefore,  the  interpolation  will  be  performed  only  in  the  horizontal 
direction  using  one  dimensional  cubic  spline  functions. 

Define  the  following  approximations  for  the  time  derivatives  and  time 
average  respectively 


dHF       F+-F 


(4.34) 


dt 


2Af 


~       F+  +  F 
F  =  


(4.35) 


Use  of  approximations  (4.34)  and  (4.35)  into  (4.28)-(4 .31 )  gives 


U    -U 

2Af 


rx  -  0, 


(4.36) 


v   -v        1 

+  — 


dip*      d<p 


2At       2 


{  fy      fy 


+  r,-0, 


(437) 
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^\^1   .0,  (4.38a) 


dy        dz 


*L  +  !t!L  «  o,  <438b> 


dy        dz 


r-e-    i  ae*/  + 
2Ar     2    az  v 


-(w+  +  w-)  +  r,  -  0.  <4^) 


Group  the  terms  at  time  (t-At)  and  define  new  auxiliary  variables  pl ; 


Px   '   -«", 


(4.40) 


(4.41) 


p3  .  -e'  +  Ar-^w".  W.42) 

dz 


28 


<7,; 


Similarly,  group  the  terms  at  tune  ft  *■  A  /Mod  define  the  auxdiarv  variables 


ql  =  u 


(4  43) 


ay 


(4.44) 


^3   -   e^  +  Af  — W*.  (4.45) 

dz 


Use  of  (4.40)-(4.45)  allows  to  rewrite  (4.36),  (4.37)  and  (4.39)  as  foiJows: 
q{   -    -£,-2Af  r,  ,  (4.46) 


q2  =  -p2-2At  r2,  (4.47) 


?3  -  -p3-2Ar  r3.  (4  48) 


The  variables  /?,  and  rk  can  be  calculated  explicitly  and  with  their  values  the 
<7,  variables  can  be  obtained. 
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To  calculate  the  values  of  the  variables  u.  v,  w.  6  and  q>  at  the  time  level 
(t  t  At)  it  is  necessary  to  solve  the  system  of  equations  composed  by  (4.38b),  (4.43)- 
(4.45)  and  the  hydrostatic  equation  (4  24b)   The  procedure  used  by  Robert  et  al.(  1985) 
consists  of  first  solving  an  elliptic  equation  for  q> 

[n  this  model  it  is  not  possible  to  state  exact  boundary  conditions  for  <p,  due 
to  the  use  of  the  rigid  lid  assumption.  However,  the  w  field  has  the  exact  boundary 
conditions  (3.1 3).  The  follow  ing  elliptic  equation  for  w  can  be  derived  from  the  system 
of  equations  at  time  level  ft  +  At)  : 

&w*  +  A^g    ar    &mT_  m    Atg    &<h  +  fQ2_         (4  49) 
dz2         %       dz      dy2  %      dy2      dydz  ' 

Once  w  was  calculated,  8  was  determined  using  (4.45),  and  — 5_  and  v  were 

dy 
obtained  using  the  procedure  described  in  Chapter  III.  After  adjusting  the  v  field  to 

satisfy  <  v>  -  0,  a  new  -^5-  field  was  calculated  using  the  new  values  of  v  in  (4.44) 

ay 

and  the  w  field  was  also  adjusted  to  satisfy  the  continuity  equation  (4.38a). 

Finally,  both  in  the  explicit  and  the  semi-implicit  schemes  the  spatial 
derivatives  are  approximated  by  centered  in  space  finite  differences.  Following  Monk 
( 1989),  the  derivatives  are  first  evaluated  at  the  grid  points  and,  after  that  the  values  are 
interpolated  for  the  position  of  the  fluid  particle.  This  procedure  is  computationally 
more  efficient  than  the  calculation  of  the  derivatives  at  the  fluid  particle  points,  which 
would  double  the  number  of  interpolations. 
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V.   EXPERIMENTS  AND  RESULTS 


A       PARAMETERS  USED  FN  THE  EXPERIMENTS 

The  following  values  for  the  constants  were  used  in  all  experiments: 
L  =  3600  km, 
H  =  12  km. 
A=  12°  K, 

e .  =  ^oo°  k, 

g  ^  9.81  ms\ 
f-  10  V, 
DH  =  10V  , 


ae 


dz 


9  km, 

£  -  4°  K  km "' 


For  the  semi-Lagrangian  semi-implicit  (SLSI)  case  there  were  two  control  runs 
using  b.y  -  20  km  and  Az  =  167  m.  Ail  the  remaining  experiments  used  normal 
resolution,  with  tk.y-  40  km  and  Az=  333m.  The  cross-sections  of  the  initial  fields  are 
shown  in  Figs.  2-5.  There  is  a  warm  front  with  thermally  indirect  circulation  near  y- 
900  km  and  a  cold  front  with  termally  direct  circulation  near  y-  2700  km 

In  order  to  assess  the  evolution  of  the  frontogenesis  process,  the  following 
parameter  is  defined  at  the  lowest  computational  level: 

lAei 
ae  <5» 


d  - 


dy 


max 
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Figure  2.  Initial  u  component  field  (m/s). 
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Figure  3.  Initial  v component  field  (m/s). 
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Figure  4.  Initial  w  component  field  (m/s). 


Figure  5.  Initial  e  field  (°K). 
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where  A6  is  the  maximum  variation  of  0  along  the  y  direction  This  parameter,  used 
by  Williams  et  al.(1991),  gives  a  reasonable  measure  of  the  width  of  the  frontal  zone 
During  the  first  experiment  wit!  he  semi-implicit  scheme,  a  high  frequency 
oscillation  was  observed  in  the  frontal  width.  That  noise  was  not  present  in  the 
experiments  carried  out  with  the  semi-Lagrangian  explicit  (SLEX)  scheme  In  order  to 
eliminate  the  noise,  a  time  filtering  technique  similar  to  the  one  described  by  Asselin 
( 1972)  was  used.  The  time  filtering  of  a  variable  Q  is  defined  as 


Q(t)  -  Q(thvlQ(t+At)-2Q(tyQ(t-At)),  <52) 

where  Q(t+A  />has  been  previously  obtained  by  using  the  respective  predictive  equation 
with  the  unaveraged  value  of  Q(t).  The  time  filter  defined  above  is  expected  to  affect  the 
higher  frequencies  only. 

Several  experiments  were  performed  with  the  objective  of  assessing  how  the  use 
of  time  filtering  would  affect  the  solutions,  as  will  be  described  next. 

B.      DESCRIPTION  OF  THE  EXPERIMENTS 

The  experiments  basically  consisted  of  carrying  out  time  integrations  of  the  model 
using  the  SLEX  and  SLSI  schemes,  starting  from  the  same  initial  conditions  and  using 
different  combinations  of  time  step  A /and  filtering  factor  y.  The  experiments  were 
intended  to  give  information  about  the  following  points: 

•  How  well  the  different  schemes  handle  the  formation  of  large  gradients; 

•  How  the  use  of  time  filtering  affects  the  final  solutions;  and 

•  How  the  solutions  are  affected  by  the  use  of  longer  time  steps. 


34 


Table  1  shows  the  experiments  performed  with  the  SLEX  scheme,  as  a  function 
of  time  step  A/1  and  filtering  factor  y  The  experiments  carried  out  using  the  SLSI 
scheme  are  identified  on  Table  2. 

TABLE  1    EXPERIMENTS  WITH  THE  SLEX  SCHEME 




A  t  (sec.) 

y  =0.0 

y  =  0.01 

Y  -  0.03 

180 

SLEX01 

SLEX02 

SLEX03 

360 

SLEX04 

SLEX05 

SLEX06 

600 

SLEX07 



SLEX08 

TABLE  2.   EXPERIMENTS  WITH  THE  SLSI  SCHEME. 


1  At  (sec.) 

Y  =0.0 

Y=0.01 

y=0.03 

y=0.05 

Y=0.07 

Y=0.1 

1    i8° 

SLSI01 

SLSI02 

SLSI03 

SLSI23 

— 

360 

SLSI04 

SLSI05 

SLSI06 

SLSI24 

— 

J   600 

SLSI07 

SLSI08 

SLSI09 

SLSI25 

— 

1200 

— 

SLSI  10 

SLSI11 

SLSI  12 

— 

1800 

— 

SLSI13 

SLSI14 

SLSI  15 

— 

2400 

— 

SLSI  16 

SLSI17 

SLSI18 

SLSI  19 

3600 

— 

SLSI20 

SLSI21 

SLSI22 

SLSI26 

SLSI27 
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There  were  two  high-resolution  experiments,  both  using  the  SLSI  scheme,  with 
tune  steps  of  180  sees  In  the  first  of  them  ,  identified  by  SLHR01,  the  time  filtering 
factor  y  =  0.0 1  was  used  .  and  in  the  second  one  (SLHR02)  y  was  set  equal  to  0.05.  The 
velocity  components  and  temperature  fields  at  tune  /  -  40  hours,  for  experiment 
SLHR01  are  shown  in  Figs.  6-9  The  occurrence  of  frontogenesis  for  the  surface  cold 
front  and  frontolysis  for  the  warm  front  can  be  seen.  These  processes  are  due  to  the 
relation  between  the  wind  deformation  and  temperature  disturbance  fields,  that  gives 
convergence  in  the  cold  front  region  and  divergence  for  the  warm  front.  The  wand  v 
fields  give  a  direct  circulation  around  the  cold  frontal  zone  and  an  indirect  circulation 
about  the  warm  frontal  zone.   The  u  field  develops  cyclonic  shear  at  the  cold  front. 

The  evolution  of  the  frontal  width  (d-vaJue)  for  SLHR01  is  presented  in  Fig.  10. 
An  almost  linear  decrease  in  the  frontal  scale  can  be  seen.  This  curve  is  different  from 
that  obtained  by  Williams  et  al.(1991),  because  in  that  case  there  was  a  momentum 
diffusion  term  that  gave  an  asymptotic  behavior  in  the  evolution  of  the  frontal  width 
towards  a  balance  condition.  In  the  present  model  a  continuous  reduction  of  the  frontal 
scale  is  expected  until  a  limit,  when  the  grid  resolution  is  reached.  Figure  1 1  shows  the 
evolution  of  the  d-vaJue  for  experiment  SLEX01 .  It  can  be  seen  that  the  minimum  value 
of  the  frontal  width  is  reached  at  approximately  t  -  40  hours.  Note  that  this  curve  is 
more  linear  than  that  obtained  in  experiment  SLHR01.  Also,  the  d- value  curves 
obtained  for  experiments  SLEX02  and  SLEX03  were  virtually  coincident  with  the  one 
for  SLEX01 .  This  result  shows  that  the  use  of  stronger  filtering  effect  (larger  y )  does  not 
affect  the  evolution  of  the  frontogenesis  process  represented  by  the  SLEX  scheme. 

During  the  first  experiment  with  the  SLSI  scheme  (SLSI01)  the  d-va/ue  curve, 
shown  in  Fig.  12,  presented  a  high  frequency  oscillation.  The  time  integration  could  be 
carried  out  until  time  t  -  36  hours  and  after  that,  the  iterative  procedures  no  longer 
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Figure  6.  u  component  field  (m/s)  for  SLHR01  at  t  =  40  hours. 
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Figure  7.  v  component  field  (m/s)  for  SLHR01  at  t  =  40  hours. 
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Figure  8.  W  component  field  (m/s)  for  SLHR01  at  t  =  40  hours. 


Figure  9.  6  field  (°K)  for  SLHR01  at  t  =  40  hours. 
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Figure  10.  Evolution  of  the  frontal  width  (d-value)  for  SLHR01. 
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Figure  11.  Evolution  of  the  frontal  width  (d-value)  for  SLEX01. 
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achieved  convergence  and  the  integration  was  terminated.  It  is  important  to  note  that 
the  irregular  behavior  of  the  curve  was  not  related  to  continuous  amplification  of  the 
solution  due  to  numerical  instability  and  that  the  curve  presents  a  general  tendency  to 
the  frontogenesis  solution  obtained  in  the  SLEX  experiments.  The  high  resolution 
experiment  SLHR01  also  shows  a  low  frequency  oscillation  about  its  linear  decrease  in 
d 


i — | — i — | — i — | — i — \ — i      |      i      |      i      |      r 
0.0       0.2       0.4       0.6       0.8        1.0        1.2        1.4        1.6        1.8       2.0 

Time  (days) 


Figure  12.  Evolution  of  the  frontal  width  (d-value)  for  SLSI01. 

Ln  order  to  eliminate  the  noise  the  time  filter  described  in  the  last  section  was 
introduced.  Figure  13  shows  the  evolution  of  the  frontal  width  obtained  in  experiment 
SLSI02,  with  the  filtering  factor  y  set  equal  to  0.01 .  It  can  be  seen  from  the  figure  that 
although  the  filtering  effect  was  small,  it  was  adequate  to  suppress  the  high  frequency 
oscillations. 
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Figure  13.  Evolution  of  the  frontal  width  (d-value)  for  experiment  SLSI02. 

Figures  1 1  and  13  show  that  both  for  the  the  SLSI  and  SLEX  experiments  the 
minimum  frontal  width  is  reached  approximately  at  time  t-  40  hours  and  after  that  the 
d-value  curves  oscillate.  Those  oscillations  can  be  related  to  the  fact  that  the  front 
reached  a  width  corresponding  to  the  model  resolution  and  after  that  point  the  solutions 
were  no  longer  physically  consistent. 

Tables  3  and  4  show  the  minimum  value  of  the  frontal  width  for  the  SLEX 
experiments  with  time  steps  180,  360  and  600  sec.,  and  the  corresponding  SLSI 
experiments  for  which  the  time  integrations  could  be  carried  out  until  t  =  40  hours. 
Also  presented  is  the  time  when  the  minimum  d-value  ocurred. 

Comparison  of  the  values  presented  in  Tables  3  and  4  show  that  the  minima  d- 
values  obtained  in  the  SLSI  experiments  are  smaller  than  those  obtained  for  the 
corresponding  SLEX  experiments  with  same  A  /and  y .  Also,  in  the  SLEX  experiments 
the  minimum  frontal  width  increases  as  A  /increase  whereas  in  the  SLSI  case  there  is  no 
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significant  change  in  the  minimum  d-value  due  to  changes  in  the  time  step  It  is 
interesting  to  note  that  for  the  SLEX  case  a  shorter  tune  integration  is  necessary  for  the 
minimum  to  occur  for  increasing  tune  step,  while  the  inverse  tendency  is  observed  in  the 
SLSI  case  that  is.  a  longer  time  integration  is  necessary  for  the  minimum  frontal  width 
to  occur,  for  larger  Ar 

TABLE  3.  MINIMUM  D-VALUE  FOR  SLEX  EXPERIMENTS. 


A  t  (sec.) 

Y 

Experiment 

Minimum 
d-value  (km) 

Time  (hour) 

180 

0.0 

SLEX01 

49.0 

455 

180 

0.01 

SLEX02 

49.0 

45.5 

180 

0.03 

SLEX03 

49.0 

45.5 

360 

0.0 

SLEX04 

68.4 

39.6 

360 

0.01 

SLEX05 

68.3 

39.9 

360 

0.03 

SLEX06 

683 

396 

600 

0.0 

SLEX07 

69.9 

39, 

600 

0.03 

SLEX08 

69.8 

39.5 

Therefore,  it  can  be  concluded  that  the  SLSI  scheme  gives  a  better  representation 
of  the  scale  collapse  process  than  that  obtained  from  the  SLEX  scheme,  since  narrower 
frontal  widths  can  be  resolved  by  the  former  than  by  the  latter. 
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TABLE  4.   MINIMUM  D-VALUH  FOR  SLSI  EXPERIMENTS 


A  t  (sec.) 

Y 

Experiment 

Muimum 
d-value  (km) 

Time 

(hours) 

180 

0.01 

SLSI02 

43.4 

41.5 

180 

0.03 

SLSI03 

43.4 

41.7 

360 

001 

SLSI05 

43.4 

41  5 

360 

0.03 

SLSI06 

43.4 

42.0 

600 

0.03 

SLSI09 

43.1 

42.0 

Experiment  SLSI03  used  y=  0.03  and  the  d-vaJuecuvxz  obtained  was  practically 
coincident  with  the  one  obtained  in  SLSI02,  with  y=0.01. 

In  experiments  SLSI04  and  SLSI07  no  filtering  was  used  (y=0.0)  and 
approximately  at  time  t  -  37  hours  the  numerical  integration  was  interrupted  because 
the  solutions  were  no  longer  able  to  achieve  convergence  in  the  iterative  procedures 
However,  the  use  of  larger  values  of  y  (SLSI05.  SLSI06,  SLSI24,  SLSI08.  SLSI09. 
SLSI25)  allowed  the  normal  execution  of  the  time  integrations.  Figure  14  shows  the  d- 
vaJue  curves  for  experiments  SLSI08  and  SLSI09.  It  can  be  seen  that  the  oscillations 
were  not  completely  filtered  out  with  y =0.01  and  a  stronger  filtering  effect  (y =0.03)  was 
required.  The  experiments  with  longer  time  steps  showed  that  as  the  time  step  increased 
a  larger  value  for  y  was  required  for  an  appropriate  filtering  of  the  noise.  Also,  all 
experiments  with  the  SLSI  scheme  showed  that  the  largest  amplitudes  of  the  noise 
occurred  during  the  first  24  hours  of  integration.    Furthermore,  those  experiments 
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showed  that  the  noise  was  not  related  to  imbalance  in  the  initial  conditions,  since  the 
mass  and  wind  fields  were  initially  adjusted  using  the  quasi-geostrophic  circulation 
equation  (4.8)  and  the  oscillations  were  not  present  before  t  -  02  hours.  The  presence 
of  the  noise  was  investigated  by  examining  the  u,  v  and  6  fields  obtained  at  the  lowest 
computational  level  of  experiment  SLSI07.  At  time  t  =  13  hours,  that  corresponds 
approximately  to  maximum  amplitude  of  the  noise,  there  was  no  evidence  of  small  scale 
spatial  oscillations  in  those  fields.  Although  the  origin  of  the  noise  could  not  be 
determined,  the  experiments  showed  that  the  use  of  the  time  filter  was  effective  in 
eliminating  the  oscillations  with  relatively  small  values  of  y.  It  is  expected  that  this 
filtering  would  not  significantly  affect  the  lower  frequency  solutions. 
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Figure  14.  Evolution  of  the  fronal  width  (d-value)  for  experiments  SLSI08  and  SLSI09. 

Another  point  of  interest  in  this  study  was  the  computational  effort  employed  in 
each  of  the  experiments,  since  one  of  the  objectives  of  the  SLSI  scheme  is  to  allow  stable 
solutions  with  long  time  steps.     In  order  to  evaluate  the  relative  computational 
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efficiency,  the  CPU  tune  spent  tor  running  each  experiment  for  42  hours  was  measured 
and  normalized  by  the  CPU  tune  used  in  SLSI01.  It  was  observed  that  the  CPU  tune 
was  not  affected  by  changes  in  the  value  of  y.  so  that  the  results  are  presented  in  I  able 
5  as  a  function  of  numerical  scheme  and  tune  step  only  As  expected,  the  computational 
effort  decreases  almost  linearly  for  increasing  tune  step  Also,  for  a  given  tune  step,  the 
SLEX  scheme  is  more  efficient  than  the  SLSI  one.  This  is  observed  because  the  SLSI 
technique  requires  the  solution  of  an  elliptic  equation  every  tune  step.  It  should  also  be 
pointed  out  that  in  this  study  the  SLSI  formulation  makes  the  spatial  interpolations 
necessary  only  in  the  horizontal  direction,  whereas  in  the  SLEX  scheme  the 
interpolations  are  performed  using  bicubic  splines.  Thus,  a  formulation  of  the  SLSI 
scheme  using  two-dimensional  interpolations  could  give  a  further  increase  in  the 
computational  effort.  The  advantage  of  the  SLSI  scheme  appears  clearly  when 
integrations  with  time  steps  as  long  as  3600  sec  are  performed  with  perfectly  stable 
solutions,  while  the  SLEX  scheme  did  not  allow  time  steps  longer  than  600  seconds 
In  order  to  assess  the  accuracy  of  the  several  solutions,  a  second  high  resolution 
(SLHR02)  control  run  was  carried  out  with  the  SLSI  scheme.  The  parameters  used 
were  A/  -  180  sec.  and  y  =  0.05.  This  value  for  y  was  chosen  to  guarantee  that  the 
solutions  would  not  be  affected  by  noise.  The  solutions  of  the  control  run  were  linearly 
interpolated  for  the  normal  resolution  grid  point  positions  when  necessary,  because  of 
the  staggering  of  the  variables.  The  differences  between  the  values  obtained  for  the 
control  run  and  those  obtained  for  the  experiments  with  normal  resolution  were 
calculated  at  each  grid  point.  The  results  were  used  to  calculate  the  domain  RMS 
differences  for  each  variable.  Two  sets  of  comparisons  were  performed:  In  the  first  one 
the  velocity  components  and  temperature  fields  for  selected  SLSI  experiments  at  t-  36 
hours  were  compared  with  the  ones  obtained  from  SLHR02  at  the  same  tune    That 
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tune  was  considered  appropriate  to  give  a  well  characterized  atmospheric  front  and 
sufficiently  away  from  the  inconsistent  solutions  obtained  after  the  model  resolution 
was  achieved 

TABLES.  NORMALIZED  CPUTIME  FOR  42-HOUR  MODEL  INTEGRATION 


|  A  t  (sec.) 

SLEX 

SLSI 

180 

0.84 

1.00 

360 

0.42 

0.50 

600 

0.25 

0.30 

1200 

— 

0  16 

1800 

— 

0.11 

2400 

— 

0.08 

3600 

— 

0.06 

The  experiments  chosen  for  this  first  set  of  comparisons  were  those  with  y=  0.05 
for  all  time  steps,  except  A/=  3600  sec,  where  y=  0.07  was  used.  The  accuracy  of  the 
solutions  was  expressed  in  terms  of  RMS  differences  in  the  values  the  of  u,  v.  h  and  6, 
given  in  Table  6. 

The  values  obtained  show  that  the  u  and  6  fields  have  small  errors  as  compared 
to  the  range  of  the  values  observed  in  the  domain  (-35  to  +22  Tit  for  6  and  -50  to  + 1 5 
m  s  for  u).  On  the  other  hand,  the  errors  obtained  in  the  v  and  w  fields  were  relatively 
large  compared  with  the  magnitude  of  the  reference  values  obtained  for  SLHR02.  The 
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large  differences  can  be  associated  with  the  tact  that  the  v  and  w  components  are 
closely  related  to  the  divergence  in  the  model  and  their  values  can  be  significantly 
affected  by  the  presence  of  gravity  waves  Time  series  of  v  and  \\  were  examined  and 
they  confirmed  the  presence  o{ internal  gravity  waves  in  the  solutions  Thus,  a  single 
sample  would  not  be  sufficiently  representative  of  the  actual  level  of  accuracy  obtained 
m  those  fields. 

TABLE  6    RMS  DIFFERENCES  FOR  TIME  T=  36  HOURS 


Experiment 

At  (sec.) 

Y 

u(m/s) 

v(m/s) 

w(m/s) 
*103 

e(°K) 

SLSI23 

180 

0.05 

0.170 

0.134 

0.58 

0.087 

SLSI24 

360 

0.05 

0.180 

0.168 

0.92 

0.088 

SLSI25 

600 

0.05 

0.220 

0.256 

1.62 

0.098 

SLSI12 

1200 

0.05 

0.373 

0.526 

3.40 

0.139 

SLSI15 

1800 

0.05 

0.554 

0.796 

5.87 

0.191 

SLSI18 

2400 

0.05 

0.764 

1.038 

7.23 

0.254 

SLSI22 

3600 

0.07 

1.204 

1.270 

8.56 

0.377 

The  differences  in  the  u  and  6  fields  obtained  at  each  grid  point  for  experiments 
SLSI23  and  SLSI15  are  shown  in  Figs.  15-18  It  can  be  seen  that  the  largest  differences 
occur  close  to  the  frontal  region  for  the  u  field,  but  ney  are  almost  uniformily 
distributed  at  the  lowest  levels  in  the  0  field  for  experiment  SLSI15 
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Figure  15.  Differences  in  u  component  field  (m/s)  for  experiment  SLSI23  at  t  =  36 
hours. 


Figure  16.  Differences  in  u  component  field  (m/s)  for  experiment  SLSI15  at  t  -  36 
hours. 
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Figure  17.  Differences  in  6  field  (<>K)  for  experiment  SLSI23  at  t  =  36  hours. 


Figure  18  Differences  in  6  field  (°K)  for  experiment  SLSI15  at  t  =  36  hours 
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The  second  set  of  comparisons  used  as  reference  the  same  solutions  obtained  for 
experiment  SLHR02  at  time  t  -  36  hours.  The  frontal  width  at  that  time  in  the  high 
resolution  experiment  was  99.5  km  For  those  SLSI  experiments  where  the  time  filter 
controlled  the  high  frequency  oscillations,  the  frontal  width  d-  100  km  was  uniquely 
related  to  a  certain  time  If  the  solutions  had  no  error  we  could  suppose  that  for  a 
certain  d-value  the  velocity  components  and  temperature  fields  should  be  the  same  as 
those  obtained  from  the  control  run.  Thus,  the  geophysical  fields  for  the  tunes 
corresponding  to  a  d- valve  approximately  equal  to  100  km  were  chosen  for  the  second 
set  of  comparisons.  Table  7  gives  the  RMS  differences  obtained  and  the  corresponding 
frontal  widths  and  times  used  for  each  comparison  The  values  obtained  show  that  the 
frontogenesis  process  is  slower  for  longer  time  steps,  since  the  frontal  width  of  100  km 
is  reached  at  a  later  time  for  larger  A/.  The  RMS  differences  for  u  and  6  are  again 
relatively  small  compared  to  the  magnitude  of  the  values  obtained  in  the  reference 
solution.  The  differences  for  v  and  w  are  still  relatively  large  with  respect  to  the 
reference  values,  as  expected.  Figures  19-22  show  the  grid  point  differences  of  the  u 
and  6  fields  for  experiments  SLSI23  and  SLSI  1 5. 

The  figures  show  that  the  largest  differences  occur  close  to  the  front.  This 
characteristic  of  semi-Lagrangian  schemes  in  which  the  large  errors  areconfined  around 
the  scale  collapse  region  was  observed  by  Kuo  and  Williams  (1990).  Such  a  property 
is  desirable,  since  the  regions  away  from  the  front  will  not  undergo  a  large  effect  of 
errors  generated  close  to  the  region  where  strong  gradients  are  present. 

Figures  23-29  show  the  temperature  distribution  at  the  lowest  computational  level 
(z=  167  m)  for  the  experiments  listed  in  Table  5  (solid  line)  for  comparison  with  the 
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values  obtained  from  the  control  run  (dashed  line)    It  can  be  seen  that  as  the  time  step 

increases  the  differences  increase  but  the  largest  errors  remain  close  to  the  frontal 

region. 

TABLE  7  RMS  DIFFERENCES  FOR  d  -  100  km 


RMS  DIFFERENCES 

Experiment, 

u(mys) 

v(m/s) 

w(nvs) 

eCK) 

d(km) 

Time 

At(sec .),  y 

•10  " 

(hours) 

SLSI03, 

0.667 

0.398 

2.91 

0.215 

101.0 

.8.0 

180,  0  03 

SLSI06, 

0.683 

0.422 

3.01 

0.240 

101.4 

38.2 

360,  0.03 

SLSI09, 

0.731 

0.468 

3.24 

0.280 

1008 

38.5 

600,  0.03 

SLSI12, 

0.914 

0.593 

3.78 

0.388 

101  2 

39.3 

1200,0.05 

SLSI15, 

1.031 

0.622 

3.92 

0.452 

101.5 

39.5       1 

1800,0.05 

SLSI19, 

1.155 

0.701 

4.26 

0.531 

101.5 

40.0 

|  2400,  0.07 

I  SLSI27, 

1  197 

0815 

547 

0.588 

103 

40  0 

I  3600,0.10 
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Figure  19.  Differences  in  u  component  field  (m/s)  for  experiment  SLSI03  at  time 
corresponding  to  d- 100  km. 
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Figure  20.  Differences  in  u  component  field  (m/s)  for  experiment  SLSI15  at  time 
corresponding  to  d- 100  km. 
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Figure  21.  Differences  in  8  field  (°K)  for  experiment  SLSI03  at  time  corresponding  to 
d-lOOkm. 


Figure  22.  Differences  in  6  field  (""K)  for  experiment  SLSI15  at  tune  corresponding  to 
d-lOOkm. 
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For  this  model  the  Courant  number  a  is  calculated  using  the  following  expression: 


"  -  dvU^J    +lcJ    ) 


Af 


max 


max 


(Ay/2) 


(5.3) 


where  cymu  is  the  phase  speed  of  the  fastest  internal  gravity  wave.  An  important  result 
from  this  study  is  that  for  the  SLSI  experiments,  o  ranged  from  0.5  (Af  =  180  sec.)  to 
10.4  (A  t-  3600  sec.) ,  as  opposed  to  the  CFL  stabilty  criterion  that  would  require  o  to 
be  less  than  or  equal  to  1 ,  and  the  values  of  the  RMS  differences  obtained  for  u  and  8 
in  the  experiments  with  large  time  steps  remained  relatively  small.  This  result  indicates 
that  the  gain  in  computational  efficiency  does  not  affect  significantly  the  accuracy  of  the 
solutions  obtained  with  the  SLSI  scheme. 
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Figure  23.  Temperature  disturbance  at  z  =  167  m,  for  experiment  SLSI03  (solid  line) 
and  control  run  (dashed  line),  at  time  corresponding  to  d-100  km. 
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Figure  24.  Same  as  Fig.  23  except  for  experiment  SLSI06. 


Figure  25.  Same  as  Fig.  23  except  for  experiment  SLSI09. 


55 


-10-i 


-40- 


v* 


Control    Run 
SLSI12 


TT 


0 


I  I  I  I  I  I 
600 


M  '  ' 
1200 


i  i  i  I  i  i  I  I  i  I  i  l  I  I  I  |  I  i  i  I  I  I 
1800     2400     3000     3600 


y(  km) 


Figure  26.  Same  as  fig.  23  except  for  experiment  SLSI12. 


-i0q 

ST         - 
w  -20- 

S* 

s 

[i 

<o 

^>^ 

Y 

«-> 

'^^ 

U 

4, 

-^>^ 

X! 

1     V 

"   -30- 

1     \ 

V    \ 



--  Control    Run 

\    \ 

ci  qi  1  *i 

T- 1 1 1 

jLjI  I  D 

-40 
{ 

I   1   I   1   I   1  I   1  I   I 
D                 600 

I     I     |     I     I     I     I     I     |     l     I     l 

1 200             1 800 

y(  km) 

i         III 

2400 

riii 

3000 

1      1       1      | 

3600 

Figure  27.  Same  as  Fig.  23  except  for  experiment  SLSI15. 
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Figure  28.  Same  as  Fig.  23  except  for  experiment  SLSI19. 


-10- 


« 


E- 


-20- 


-30- 


Control    Run 
SLSI27 


40    |    i    i    i    i    i    |    i    i    i    i    i    |    i   i    i    i    i    |    i    i    i    i    i    |    i    i    i    i    i    |    i   i    i    i    i    | 
0  600  1200  1800  2400  3000  3600 

y(km) 


Figure  29.  Same  as  Fig.  23  except  for  experiment  SLSI27 
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VI.   SUMMARY  AND  CONCLUSIONS 

In  this  study  a  numerical  model  based  on  the  hydrostatic  Boussmesq  equations 
was  used  to  simulate  atmospheric  frontogenesis  driven  by  an  irrotational  non-divergent 
deformation  wind  field.  The  semi-Lagrangian  technique  was  employed  to  integrate 
numerically  the  prognostic  equations.  The  model  assumes  periodic  boundary 
conditions  in  the  cross-front  direction  and  the  rigid  lid  assumption  is  used  as  the  upper 
boundary  condition.  The  model  also  neglects  along  front  variations.  A  basic  sinusoidal 
deformation  flow  is  introduced  as  the  forcing  of  the  frontogenesis  process.  Experiments 
were  performed  using  the  semi-Lagrangian  technique  associated  with  two  different  time 
schemes:  explicit  and  semi-implicit.  In  the  semi-Lagrangian  explicit  case  (SLEX) 
bicubic  splines  were  used  to  interpolate  the  variables  in  space.  In  the  semi-Lagrangian 
semi-implicit  case  (SLSI)  the  variables  were  interpolated  in  the  y-direction  only,  using 
one-dimensional  splines.  A  frontal  width  parameter  ( d-vaJue )  was  defined  at  the  lowest 
computational  level  to  describe  the  evolution  of  the  frontogenetical  processes.  The 
experiments  with  the  SLEX  technique  were  successful  in  replicating  the  formation  of  a 
realistic  front  in  approximately  40  hours.  Different  time  steps  were  used  for  the 
integration  of  the  model  and  solutions  which  were  both  numerically  and  physically 
consistent  were  obtained  with  the  SLEX  scheme  for  values  of  Courant  numbers  as  large 
as  1.7. 

The  experiments  with  SLSI  scheme  presented  high  frequency  oscillations  that  were 
eliminated  by  using  a  time  filter.  However,  no  small  scale  spatial  oscillations  were 
observed  in  the  solutions.  Different  filtering  effects  were  tested  by  changing  the  value 
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of  the  filtering  parameter  y  The  solutions  obtained  with  the  SLEX  scheme  were  not 
affected  by  the  time  filter,  whereas  the  S LSI  solutions  showed  that  stronger  filtering  was 
necessary  for  larger  time  steps 

The  SLSI  scheme  showed  to  be  more  successful  in  handling  the  scale  collapse 
process  than  the  SLEX  scheme,  since  the  fronts  reproduced  by  the  former  had 
minimum  widths  smaller  than  those  obtained  by  the  latter  The  SLSI  scheme  presented 
a  tendency  of  slowing  down  the  frontogenesis  process  for  increasing  tune  steps 

Experiments  were  performed  with  the  SLSI  scheme  with  tune  steps  as  long  as  3600 
seconds,  corresponding  to  a  Courant  number  greater  than  10.  The  solutions  obtained 
were  perfectly  stable,  and  the  accuracy  was  not  significantly  degraded,  even  for  longer 
time  steps 

The  SLSI  solutions  also  had  the  characteristic  of  constraining  the  largest  errors 
close  to  the  scale  collapse  region.  Such  a  characteristic  is  desirable  since  the  errors  will 
have  a  smaller  impact  on  the  solutions  in  regions  away  from  the  front. 

The  results  obtained  in  this  study  suggest  that  the  SLSI  technique  is  appropriate 
for  mesoscale  regional  models  since  it  is  computationally  efficient  and  produces  accurate 
results.  A  different  formulation  of  the  scheme  where  two-dimensional  interpolation  of 
variables  were  allowed  should  be  studied  to  verify  if  there  would  be  a  significant  positive 
impact  on  the  accuracy  of  the  solutions.  Also,  effects  of  surface  topography  and 
physical  processes  like  advection  of  the  frontal  system,  friction,  vertical  shear  of  the 
basic  wind  field  and  moisture  should  be  investigated  in  future  studies 
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